R Notebook

Greenness Exposure Assessment in Environmental Health Research

Amman, July 5th 2025

Leire Luque García

University of California San Diego

EXERCISE 1: Mitigation pathway

In this exercise, we explore how different greenness metrics influence health through the mitigation pathway. As has already been noted, the mitigation pathway influences human health by reducing environmental stressors, including air pollution, heat, and noise, through the presence of surrounding vegetation.

DESCRIPTION

We conducted a prospective study involving 6,243 adult participants monitored during the warm season (May to September), to investigate the potential influence of mean daytime temperature on heart rate, and whether this association is modified by surrounding vegetation. Each participant was equipped with a validated wearable device capable of continuously recording daytime heart rates. NDVI and tree cover (%) were estimated at the residential level using remote sensing data aggregated within a 500-meter buffer surrounding each participant’s home address. In this analysis, NDVI and tree cover are treated as potential effect modifiers of the association between ambient daytime temperature and physiological stress, as reflected by heart rate.

install.packages("readr")
install.packages("dplyr")
install.packages("data.table")
install.packages("ggplot2")
install.packages("skimr")
install.packages("patchwork")
install.packages("broom")
install.packages("ggeffects")

We will start installing the required packages.

library(readr)
Mitigation <- read_csv("Mitigation.csv")

Next, load the synthetic dataset titled “Mitigation.csv”.

str(Mitigation)

Examine the structure of the dataset. The “Mitigation” dataset simulates data for 6,243 adult participants with the following variables:

• id: Indentification number
• mean_temp: Average daytime temperature (°C)
• ndvi: Normalized Difference Vegetation Index (NDVI) within a 500-m buffer of the residence
• tree_cover: Percentage of tree cover within a 500-m buffer of the residence
• mean_hr: Average daytime heart rate (beats per minute)
• sex: Sex (Female or Male)
• age: Age (years)
• ses: Socioeconomic status (I (highest) to V (lowest))

Mitigation$sex <- factor(Mitigation$sex, labels = c("Female", "Male"))
Mitigation$ses <- factor(Mitigation$ses, labels = c("I", "II", "III", "IV", "V"))

Recode categorical variables in the “Mitigation” dataset as factor.

library(skimr)

Obtain descriptive statistics for the exposure and outcome variables in the “Mitigation” dataset using “skimr”.

skim(Mitigation[, c("mean_temp", "ndvi", "tree_cover", "mean_hr")])

Check the distribution of the exposure and outcome variables.

Mitigation$ndvi_tertile <- cut(Mitigation$ndvi,
                           breaks = quantile(Mitigation$ndvi, probs = c(0, 1/3, 2/3, 1), na.rm = TRUE),
                           include.lowest = TRUE,
                           labels = c("Low", "Medium", "High"))

Mitigation$tc_tertile <- cut(Mitigation$tree_cover,
                         breaks = quantile(Mitigation$tree_cover, probs = c(0, 1/3, 2/3, 1), na.rm = TRUE),
                         include.lowest = TRUE,
                         labels = c("Low", "Medium", "High"))

Create tertiles for NDVI and tree cover.

Mitigation$ndvi_tertile <- factor(Mitigation$ndvi_tertile, levels = c("Low", "Medium", "High"))
Mitigation$tc_tertile <- factor(Mitigation$tc_tertile, levels = c("Low", "Medium", "High"))

Recode as factors.

model_mitigation <- lm(mean_hr ~ mean_temp + sex + age + ses, data = Mitigation)

This linear regression model assesses the association between mean daytime temperture (ºC) and mean heart rate (bpm), adjusting for sex, age, and socioeconomic status (SES).

summary(model_mitigation)

Interpretation
A 1°C increase in mean daytime temperature is associated with a 0.96 bpm increase in mean heart rate, adjusting for sex, age, and socioeconomic status (SES). For example, a 5°C increase would correspond to an approximate 4.8 bpm increase in heart rate. The model explains approximately 6.8% of the variance in heart rate. The residual standard error is 11.56, suggesting that predicted heart rate values deviate on average by around ±12 bpm from observed values. Overall, the model is statistically significant (p < 0.001).

Now we will assess effect modification across tertiles.

library(ggeffects)

Load the “ggeffects” package to extract predicted (marginal) effects from regression models, including models with interaction terms. These predictions are adjusted for any confounders included in the model.

model_ndvi_interaction <- lm(mean_hr ~ mean_temp * ndvi_tertile + sex + age + ses, data = Mitigation)
model_tc_interaction <- lm(mean_hr ~ mean_temp * tc_tertile + sex + age + ses, data = Mitigation)

Fit linear regression models with an interaction term.

pred_ndvi <- ggpredict(model_ndvi_interaction, terms = c("mean_temp", "ndvi_tertile"))
pred_tc <- ggpredict(model_tc_interaction, terms = c("mean_temp", "tc_tertile"))

This function calculates model-predicted heart rates across a range of temperature values, separately for each tertile.

library(patchwork)

plot_ndvi <- plot(pred_ndvi) +
  labs(
    title = "Effect Modification by NDVI",
    x = "Mean Daytime Temperature (°C)",
    y = "Predicted Mean Heart Rate (bpm)",
    color = "NDVI"
  ) +
  theme_minimal()

plot_tc <- plot(pred_tc) +
  labs(
    title = "Effect Modification by Tree Cover",
    x = "Mean Daytime Temperature (°C)",
    y = "Predicted Mean Heart Rate (bpm)",
    color = "Tree Cover"
  ) +
  theme_minimal()

plot_ndvi + plot_tc

Interpretation
The plots illustrate that the association between mean daytime temperature and heart rate is modified by surrounding vegetation. In both the NDVI and tree cover models, higher vegetation levels are associated with a progressively weaker temperature–heart rate relationship across tertiles. This pattern suggests that increased vegetation may attenuate the physiological response to heat exposure.

However, such attenuation is more pronounced in the case of tree cover. In areas with high tree cover, the slope is almost flat, indicating substantial mitigation of heat effects. In contrast, the mitigation observed for high NDVI is more modest. These results suggest that, although both metrics provide some protection, tree cover is more effective in mitigating the physiological impacts of heat.

EXERCISE 2: Restoration pathway

This exercise examines how different greenness metrics influence health through the restoration pathway. The restoration pathway describes how exposure to greenness improves health through two major theories: the Attention Restoration Theory (ART) and the Stress Reduction Theory (SRT). The SRT and ART are interconnected, grounded in biophilic principles, emphasizing the innate human connection to nature. According to the biophilia hypothesis, humans have an evolutionary tendency to seek connections with other forms of life (Wilson, 1986).

DESCRIPTION

In this exercise, you will analyze generated data from 5,000 participants to examine how greenness exposure relates to physiological stress, measured through cortisol levels. Greenness is assessed using NDVI within a 500-meter buffer around the home, while weekly time spent in green spaces (minutes) is objectively measured using individual tracking devices.

Restoration <- read_csv("Restoration.csv")

Rows: 5000 Columns: 7── Column specification ──────────────────────────────────────────────────────────────────────────
Delimiter: ","
chr (2): sex, ses
dbl (5): id, ndvi, green_minutes, cortisol, age
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

Load the synthetic dataset titled “Restoration.csv”.

str(Restoration)

spc_tbl_ [5,000 × 7] (S3: spec_tbl_df/tbl_df/tbl/data.frame)
 $ id           : num [1:5000] 1 2 3 4 5 6 7 8 9 10 ...
 $ ndvi         : num [1:5000] 0.3981 0.1379 0.1031 0.0156 0.1645 ...
 $ green_minutes: num [1:5000] 42 65 69 0 46 82 46 63 64 57 ...
 $ cortisol     : num [1:5000] 93.7 70.6 85.8 109 118.8 ...
 $ sex          : chr [1:5000] "Female" "Female" "Male" "Male" ...
 $ age          : num [1:5000] 51 48 63 42 30 47 34 36 58 43 ...
 $ ses          : chr [1:5000] "II" "III" "V" "V" ...
 - attr(*, "spec")=
  .. cols(
  ..   id = col_double(),
  ..   ndvi = col_double(),
  ..   green_minutes = col_double(),
  ..   cortisol = col_double(),
  ..   sex = col_character(),
  ..   age = col_double(),
  ..   ses = col_character()
  .. )
 - attr(*, "problems")=<externalptr> 

Examine the structure of the dataset. The “Restoration” dataset simulates data for 5,000 adult participants with the following variables:

• id: Indentification number
• ndvi: Normalized Difference Vegetation Index (NDVI) within a 500-m buffer of the residence
• green_minutes: Mean weekly time spent in green spaces (minutes)
• cortisol: Hair cortisol measures (pg/mg)
• sex: Sex (Female or Male)
• age: Age (years)
• ses: Socioeconomic status, (I (highest) to V (lowest))

Restoration$sex <- factor(Restoration$sex, labels = c("Female", "Male"))
Restoration$ses <- factor(Restoration$ses, labels = c("I", "II", "III", "IV", "V"))

Recode categorical variables in the “Restoration” dataset as factor.

skim(Restoration[, c("ndvi", "green_minutes", "cortisol")])

── Data Summary ────────────────────────
                           Values                      
Name                       Restoration[, c("ndvi", "...
Number of rows             5000                        
Number of columns          3                           
_______________________                                
Column type frequency:                                 
  numeric                  3                           
________________________                               
Group variables            None                        

Check the distribution of the exposure and outcome variables.

model_ndvi <- lm(cortisol ~ ndvi + sex + age + ses, data = Restoration)

Fit a linear regression model to examine the relationship between NDVI (exposure) and hair cortisol levels (outcome), adjusting for potential confounders (sex, age, ses).

summary(model_ndvi)

Call:
lm(formula = cortisol ~ ndvi + sex + age + ses, data = Restoration)

Residuals:
    Min      1Q  Median      3Q     Max 
-91.414 -17.447  -0.172  17.602  75.980 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept)  81.53281    1.97741  41.232  < 2e-16 ***
ndvi        -23.48901    2.74435  -8.559  < 2e-16 ***
sexMale      -4.04010    0.71633  -5.640 1.79e-08 ***
age           0.29469    0.03568   8.260  < 2e-16 ***
sesII         2.73761    1.11360   2.458    0.014 *  
sesIII        5.63268    1.12996   4.985 6.41e-07 ***
sesIV         8.92504    1.12910   7.905 3.28e-15 ***
sesV         12.86420    1.12337  11.451  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 25.29 on 4992 degrees of freedom
Multiple R-squared:  0.06367,   Adjusted R-squared:  0.06235 
F-statistic: 48.49 on 7 and 4992 DF,  p-value: < 2.2e-16

Interpretation
A 1-unit increase in NDVI (from 0 to 1) is associated with a 23.5-unit decrease in hair cortisol, adjusting for other variables. This is difficult to interpret. Using the IQR instead, allows to interpret the effect of NDVI in terms of a typical, realistic change in exposure, making the results clearer.

iqr_ndvi <- IQR(Restoration$ndvi)
Restoration$ndvi_iqr <- Restoration$ndvi/iqr_ndvi
model_ndvi_iqr <- lm(cortisol ~ ndvi_iqr + sex + age + ses, data = Restoration)

Standarize NDVI dividing it by its IQR to make the effect size more interpretable, then model its association with hair cortisol, adjusting for key covariates.

summary(model_ndvi_iqr)

Call:
lm(formula = cortisol ~ ndvi_iqr + sex + age + ses, data = Restoration)

Residuals:
    Min      1Q  Median      3Q     Max 
-91.414 -17.447  -0.172  17.602  75.980 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 81.53281    1.97741  41.232  < 2e-16 ***
ndvi_iqr    -4.30723    0.50324  -8.559  < 2e-16 ***
sexMale     -4.04010    0.71633  -5.640 1.79e-08 ***
age          0.29469    0.03568   8.260  < 2e-16 ***
sesII        2.73761    1.11360   2.458    0.014 *  
sesIII       5.63268    1.12996   4.985 6.41e-07 ***
sesIV        8.92504    1.12910   7.905 3.28e-15 ***
sesV        12.86420    1.12337  11.451  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 25.29 on 4992 degrees of freedom
Multiple R-squared:  0.06367,   Adjusted R-squared:  0.06235 
F-statistic: 48.49 on 7 and 4992 DF,  p-value: < 2.2e-16

Interpretation
A 1 IQR increase in NDVI is associated with a 4.31 pg/mg decrease in cortisol, adjusting for sex, age, and SES. The model explains 6.4% of the variance in hair cortisol levels. The residual standard error of 25.29 is relatively high, meaning that predictions from the model often deviate substantially from actual hair cortisol values. Overall, the model is statistically significant (p < 0.001).

model_time <- lm(cortisol ~ green_minutes + ses + age + sex, data = Restoration)

Fit a linear regression model to examine the relationship between hair cortisol levels (outcome) and mean weekly time spent in green spaces (exposure), adjusting for potential confounders.

summary(model_time)

Call:
lm(formula = cortisol ~ green_minutes + ses + age + sex, data = Restoration)

Residuals:
    Min      1Q  Median      3Q     Max 
-93.191 -17.135  -0.321  17.252  81.090 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)   86.31237    2.07579  41.580  < 2e-16 ***
green_minutes -0.14979    0.01335 -11.223  < 2e-16 ***
sesII          1.91357    1.10992   1.724 0.084758 .  
sesIII         4.19615    1.13207   3.707 0.000212 ***
sesIV          6.88432    1.13954   6.041 1.64e-09 ***
sesV          10.05635    1.14318   8.797  < 2e-16 ***
age            0.27694    0.03556   7.787 8.26e-15 ***
sexMale       -3.93888    0.71249  -5.528 3.40e-08 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 25.16 on 4992 degrees of freedom
Multiple R-squared:  0.07331,   Adjusted R-squared:  0.07201 
F-statistic: 56.41 on 7 and 4992 DF,  p-value: < 2.2e-16

Interpretation
A 1-minute increase in the mean weekly time spent in green spaces is associated with a 0.15 pg/mg decrease in hair cortisol, adjusting for SES, age, and sex. For example, an increase of 30 minutes per week would correspond to a 4.5 pg/mg reduction in cortisol levels. The model explains 7.3% of the variance in hair cortisol levels. The residual standard error is 25.16, indicating that predicted cortisol values often deviate notably from observed values. Overall, the model is statistically significant (p < 0.001).

library(patchwork)

Load the required package to combine ggplot objects.

p1 <- ggplot(Restoration, aes(x = ndvi, y = cortisol)) +
  geom_point(alpha = 0.3, color = "#145A32") +
  geom_smooth(method = "lm", color = "#27AE60", fill = "#ABEBC6") +
  labs(
    title = "NDVI vs. Hair Cortisol",
    x = "NDVI (Residential Greenness)",
    y = "Hair Cortisol (pg/mg)"
  ) +
  theme_minimal(base_size = 10)

Create a scatterplot showing the association between NDVI and hair cortisol levels.

p2 <- ggplot(Restoration, aes(x = green_minutes, y = cortisol)) +
  geom_point(alpha = 0.3, color = "#2C3E50") +
  geom_smooth(method = "lm", color = "#2980B9", fill = "#85C1E9") +
  labs(
    title = "Green Space Time vs. Hair Cortisol",
    x = "Green Space Time (minutes/week)",
    y = "Hair Cortisol (pg/mg)"
  ) +
  theme_minimal(base_size = 10)

Create a second scatterplot on the association between mean weekly time spent in green spaces and hair cortisol levels.

combined_plot <- p1 + p2
combined_plot

Combine plots for comparison using patchwork.

Interpretation
While both “NDVI” and “Green Space Time” are significantly associated with lower hair cortisol, “Green Space Time” shows a stronger statistical association and slightly better model fit. This suggests that while NDVI is a widely available and straightforward greenness metric, the restoration pathway likely involves direct interaction with nature. Therefore, studies should consider not only spatial availability of greenness but also the frequency and duration of individual exposure.

EXERCISE 3: Instoration pathway

This exercise examines how different greenness metrics influence health through the instoration pathway. The instoration or capacity building pathway, highlights how natural environments promote health-related behaviors by encouraging physical activity and social cohesion.

Green spaces often provide a secure, accessible, and attractive setting in which to be physically active. In this regard, numerous studies have found positive associations between increased exposure to green environments and higher levels of physical activity (Gianfredi et al., 2021), which, in turn, is linked to improved physical (Bull et al., 2020) and mental health (Rodriguez-Ayllon et al., 2019).

DESCRIPTION

In this exercise, you will analyze generated data from 5,000 participants to examine how proximity to green spaces relates to physical activity, measured as weekly minutes using individual tracking devices. Greenness exposure is assessed using two different distance measures based on the residential address:

• Straight-line distance to the closest green space (in meters)
• Network distance to the closest green space (in meters)

Instoration <- read_csv("Instoration.csv")

Rows: 5000 Columns: 7── Column specification ──────────────────────────────────────────────────────────────────────────
Delimiter: ","
chr (2): sex, ses
dbl (5): id, straight_distance, network_distance, weekly_pa, age
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

Load the synthetic dataset titled “Instoration.csv”.

str(Instoration)

spc_tbl_ [5,000 × 7] (S3: spec_tbl_df/tbl_df/tbl/data.frame)
 $ id               : num [1:5000] 1 2 3 4 5 6 7 8 9 10 ...
 $ straight_distance: num [1:5000] 433 62 296 341 233 214 188 266 148 503 ...
 $ network_distance : num [1:5000] 904 257 341 682 257 279 377 755 599 791 ...
 $ weekly_pa        : num [1:5000] 292 295 303 305 271 296 315 262 278 280 ...
 $ sex              : chr [1:5000] "Male" "Male" "Male" "Male" ...
 $ age              : num [1:5000] 36 50 43 58 60 38 52 68 54 51 ...
 $ ses              : chr [1:5000] "V" "I" "III" "V" ...
 - attr(*, "spec")=
  .. cols(
  ..   id = col_double(),
  ..   straight_distance = col_double(),
  ..   network_distance = col_double(),
  ..   weekly_pa = col_double(),
  ..   sex = col_character(),
  ..   age = col_double(),
  ..   ses = col_character()
  .. )
 - attr(*, "problems")=<externalptr> 

Examine the structure of the dataset. The “Instoration” dataset simulates data for 5,000 adult participants with the following variables:

• id: Indentification number
• straight_distance: Straight-line distance to the nearest green space (meters) based on the residence
• network_distance: Network distance to the nearest green space (meters) based on the residence
• weekly_pa: Weekly physical activity (minutes)
• sex: Sex (Female or Male)
• age: Age (years)
• ses: Socioeconomic status, (I (highest) to V (lowest))

Instoration$sex <- factor(Instoration$sex, labels = c("Female", "Male"))
Instoration$ses <- factor(Instoration$ses, labels = c("I", "II", "III", "IV", "V"))

Recode categorical variables in the “Instoration” dataset as factor.

skim(Instoration[, c("straight_distance", "network_distance", "weekly_pa")])

── Data Summary ────────────────────────
                           Values                      
Name                       Instoration[, c("straight...
Number of rows             5000                        
Number of columns          3                           
_______________________                                
Column type frequency:                                 
  numeric                  3                           
________________________                               
Group variables            None                        

Check the distribution of the exposure and outcome variables.

model_straight <- lm(weekly_pa ~ straight_distance + sex + age + ses, data = Instoration)

Fit a linear regression model to examine the relationship between straight-line distance to the nearest green space (exposure) and weekly physical activity time (outcome), adjusting for potential confounders (sex, age, ses).

summary(model_straight)

Call:
lm(formula = weekly_pa ~ straight_distance + sex + age + ses, 
    data = Instoration)

Residuals:
    Min      1Q  Median      3Q     Max 
-97.075 -20.316  -0.386  19.883 119.331 

Coefficients:
                    Estimate Std. Error t value Pr(>|t|)    
(Intercept)       300.585497   2.248761 133.667   <2e-16 ***
straight_distance  -0.023692   0.001754 -13.506   <2e-16 ***
sexMale            -1.581549   0.841731  -1.879   0.0603 .  
age                -0.062596   0.037653  -1.662   0.0965 .  
sesII               0.334750   1.311055   0.255   0.7985    
sesIII             -0.119385   1.312486  -0.091   0.9275    
sesIV               1.203676   1.320661   0.911   0.3621    
sesV                0.424646   1.324429   0.321   0.7485    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 29.72 on 4992 degrees of freedom
Multiple R-squared:  0.03637,   Adjusted R-squared:  0.03502 
F-statistic: 26.91 on 7 and 4992 DF,  p-value: < 2.2e-16

Interpretation
A 100-meter increase in straight-line distance to the nearest green space is associated with a 2.37-minute decrease in weekly physical activity, adjusting for sex, age, and socioeconomic status. The model explains 3.6% of the variance in weekly physical activity. The residual standard error indicates a substantial unexplained variability of around 30 minutes. Although the association is statistically significant (p < 0.001), the overall effect size is small, and the explanatory power is limited.

model_network <- lm(weekly_pa ~ network_distance + sex + age + ses, data = Instoration)

Fit a linear regression model to examine the relationship between network distance to the nearest green space (exposure) and weekly physical activity time (outcome), adjusting for potential confounders (sex, age, ses).

summary(model_network)

Call:
lm(formula = weekly_pa ~ network_distance + sex + age + ses, 
    data = Instoration)

Residuals:
    Min      1Q  Median      3Q     Max 
-97.596 -20.786  -0.414  20.064 122.487 

Coefficients:
                   Estimate Std. Error t value Pr(>|t|)    
(Intercept)      306.284376   2.384980 128.422   <2e-16 ***
network_distance  -0.021386   0.001513 -14.132   <2e-16 ***
sexMale           -1.668671   0.840360  -1.986   0.0471 *  
age               -0.063695   0.037591  -1.694   0.0903 .  
sesII              0.347520   1.308872   0.266   0.7906    
sesIII            -0.161458   1.310314  -0.123   0.9019    
sesIV              1.167341   1.318413   0.885   0.3760    
sesV               0.205750   1.322074   0.156   0.8763    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 29.67 on 4992 degrees of freedom
Multiple R-squared:  0.03958,   Adjusted R-squared:  0.03823 
F-statistic: 29.39 on 7 and 4992 DF,  p-value: < 2.2e-16

Interpretation
A 100-meter increase in network distance to the nearest green space is associated with a 2.14-minute decrease in weekly physical activity, adjusting for sex, age, and socioeconomic status. The model explains 4.0% of the variance in weekly physical activity. The residual standard error indicates a substantial unexplained variability of around 30 minutes. Although the association is statistically significant (p < 0.001), the overall effect size is small, and the explanatory power is limited.

p1 <- ggplot(Instoration, aes(x = straight_distance, y = weekly_pa)) +
  geom_point(alpha = 0.3, color = "#ff70a6") +
  geom_smooth(method = "lm", color = "#ff0054", fill = "#F796BB") +
  labs(
    title = "Straight-line vs. Physical Activity",
    x = "Straight-line Green Space (meters)",
    y = "Physical Activity (minutes/week)"
  ) +
  theme_minimal(base_size = 10)

Create a scatterplot showing the association between straight-line distance to the closest green space and weekly physical activity.

p2 <- ggplot(Instoration, aes(x = network_distance, y = weekly_pa)) +
  geom_point(alpha = 0.3, color = "#E4B0EE") +
  geom_smooth(method = "lm", color = "#BE39D8", fill = "#D398DF") +
  labs(
    title = "Network vs. Physical Activity",
    x = "Network Green Space (meters)",
    y = "Physical Activity (minutes/week)"
  ) +
  theme_minimal(base_size = 10)

Create a scatterplot showing the association between network distance to the closest green space and weekly physical activity.

combined_plot <- p1 + p2
combined_plot

Combine plots for comparison using patchwork.

Interpretation
While both “straight-line distance” and “network distance” to the nearest green space are significantly associated with lower weekly physical activity, “network distance” shows a slightly stronger statistical association (t = –14.13 vs. –13.51) and better model fit (R² = 4.0% vs. 3.6%). These results suggest that network distance offers a more realistic estimation of green space accessibility, as it is constructed based on actual routes.

EXERCISE 4: Biodiversity pathway

This exercise focuses on the biodiversity pathway. Biodiversity is a fundamental determinant of health, influencing both physical and psychological well-being. While biodiversity generally benefits human health, it also presents associated risks, including exposure to airborne allergens and volatile organic compounds that may trigger allergic reactions (Marselle et al., 2021).

DESCRIPTION

You will analyze generated data from 4,730 participants to examine how vegetation diversity is associated with allergic rhinitis symptoms (e.g., sneezing, nasal congestion, itchy eyes), measured as the total number of antihistamine pills taken per year.

Vegetation diversity is modeled as a continuous index based on the Shannon diversity index, which accounts for both abundance and evenness of the species present. The index ranges from 0 (no diversity) to approximately 2.5 (high diversity), with higher values indicating greater ecological diversity in the surrounding plant environment.

library(readr)
Biodiversity <- read_csv("Biodiversity.csv")

Rows: 4730 Columns: 7── Column specification ──────────────────────────────────────────────────────────────────────────
Delimiter: ","
chr (3): urbanicity, sex, ses
dbl (4): id, veg_diversity, allergy_med, age
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

Load the synthetic dataset titled “Biodiversity.csv”.

str(Biodiversity)

spc_tbl_ [4,730 × 7] (S3: spec_tbl_df/tbl_df/tbl/data.frame)
 $ id           : num [1:4730] 1 2 3 4 5 6 7 8 9 10 ...
 $ veg_diversity: num [1:4730] 1.91 1.37 1.9 2.07 1.56 2.23 1.33 1.46 1.95 1.92 ...
 $ urbanicity   : chr [1:4730] "Suburban" "Urban" "Urban" "Suburban" ...
 $ allergy_med  : num [1:4730] 4 7 7 8 5 12 6 5 11 7 ...
 $ sex          : chr [1:4730] "Male" "Female" "Male" "Male" ...
 $ age          : num [1:4730] 53 40 51 52 39 70 59 44 49 30 ...
 $ ses          : chr [1:4730] "I" "V" "III" "III" ...
 - attr(*, "spec")=
  .. cols(
  ..   id = col_double(),
  ..   veg_diversity = col_double(),
  ..   urbanicity = col_character(),
  ..   allergy_med = col_double(),
  ..   sex = col_character(),
  ..   age = col_double(),
  ..   ses = col_character()
  .. )
 - attr(*, "problems")=<externalptr> 

Examine the structure of the dataset. The “Biodiversity” dataset simulates data for 4,730 adult participants with the following variables:

• id: Indentification number
• veg_diversity: Vegetation diversity based on the Shannon diversity index, estimated within a 300-meter buffer around each participant’s residence.
• urbanicity: Urbanization level (Rural, Suburban or Urban) of the neighborhood area.
• allergy_med: The total number of antihistamine pills consumed over the the past year, as a proxy of allergic rhinitis symptom burden.
• sex: Sex (Female or Male)
• age: Age (years)
• ses: Socioeconomic status (I (highest) to V (lowest))

Biodiversity$urbanicity <- factor(Biodiversity$urbanicity, labels = c("Rural", "Suburban", "Urban"))
Biodiversity$sex <- factor(Biodiversity$sex, labels = c("Female", "Male"))
Biodiversity$ses <- factor(Biodiversity$ses, labels = c("I", "II", "III", "IV", "V"))

Recode categorical variables in the “Biodiversity” dataset as factor.

library(skimr)
skim(Biodiversity[, c("veg_diversity", "allergy_med")])

── Data Summary ────────────────────────
                           Values                      
Name                       Biodiversity[, c("veg_div...
Number of rows             4730                        
Number of columns          2                           
_______________________                                
Column type frequency:                                 
  numeric                  2                           
________________________                               
Group variables            None                        

Check the distribution of the exposure and outcome variables.

model_biodiversity <- lm(allergy_med ~ veg_diversity + sex + age + ses, data = Biodiversity)

Fit a linear regression model to examine the relationship between vegetation diversity (exposure) and antihistamine pills per year (outcome), adjusting for potential confounders (sex, age, ses).

summary(model_biodiversity)

Call:
lm(formula = allergy_med ~ veg_diversity + sex + age + ses, data = Biodiversity)

Residuals:
    Min      1Q  Median      3Q     Max 
-9.1043 -2.5355 -0.3814  2.2014 15.2418 

Coefficients:
                Estimate Std. Error t value Pr(>|t|)    
(Intercept)   -1.8879709  0.3396285  -5.559 2.86e-08 ***
veg_diversity  4.5195009  0.1353776  33.384  < 2e-16 ***
sexMale        0.1948230  0.0994269   1.959   0.0501 .  
age            0.0072995  0.0043829   1.665   0.0959 .  
sesII         -0.0001596  0.1551195  -0.001   0.9992    
sesIII        -0.2834895  0.1545145  -1.835   0.0666 .  
sesIV         -0.1052469  0.1560432  -0.674   0.5000    
sesV          -0.2214779  0.1584723  -1.398   0.1623    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 3.418 on 4722 degrees of freedom
Multiple R-squared:  0.193, Adjusted R-squared:  0.1918 
F-statistic: 161.3 on 7 and 4722 DF,  p-value: < 2.2e-16

Interpretation
A unit increase in vegetation diversity is associated with an estimated increase of 4.52 antihistamine pills per year, adjusting for sex, age, and socioeconomic status (p < 0.001). The model explains 19.3% of the variance in antihistamine use, with a residual standard error of 3.42 pills, which indicates a moderate unexplained variability.

library(dplyr)

Attaching package: ‘dplyr’

The following objects are masked from ‘package:data.table’:

    between, first, last

The following objects are masked from ‘package:stats’:

    filter, lag

The following objects are masked from ‘package:base’:

    intersect, setdiff, setequal, union

library(broom)

model_coefficients_urbanicity <- Biodiversity %>%
  group_by(urbanicity) %>%
  do(tidy(lm(allergy_med ~ veg_diversity + sex + age + ses, data = .))) %>%
  filter(term == "veg_diversity")  

Fit the previous linear regression model, but this time stratified by urbanization level. This approach allows to examine whether the strength or direction of the association between vegetation diversity and allergic rhinitis symptoms varies across rural, suburban, and urban residential contexts.

print(model_coefficients_urbanicity)

Interpretation
When stratified by urbanization level, the association between vegetation diversity and antihistamine use differs in magnitude. In urban areas, a unit increase in vegetation diversity is associated with an estimated 4.47 pill increase per year (p < 0.001), after adjusting for sex, age, and socioeconomic status. In suburban areas, the increase corresponds to 1.41 pills per year (p < 0.001), while in rural areas, the association is smaller at 0.56 pills per year (p < 0.001). These results suggest that the health impact of vegetation diversity may be modified by urbanization level.

ggplot(Biodiversity, aes(x = veg_diversity, y = allergy_med, color = urbanicity)) +
  geom_point(alpha = 0.4) +
  geom_smooth(method = "lm", se = FALSE, size = 1.2) +
  scale_color_manual(
    values = c("Rural" = "#0BDA51", "Suburban" = "#1E90FF", "Urban" = "#FF3800")  # you can change these
  ) +
  labs(
    title = "Vegetation Diversity vs. Allergy Medication",
    x = "Vegetation Diversity",
    y = "Allergy Medication (pills/year)",
    color = "Urbanization"
  ) +
  theme_minimal()

Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
Please use `linewidth` instead.

Interpretation
The relationship between urbanization and allergic symptoms (e.g.,asthma, allergic rhinitis, eczema, food allergies) has been widely studied, and epidemiological evidence consistently shows higher prevalence and severity of allergic diseases among individuals residing in urban areas when compared to those that reside in rural settings.This association may be partly attributed to reduced microbial exposure in urban environments. Urbanization has been associated to lower microbial diversity, which is considered an important factor in the development of allergic sensitization. According to the hygiene hypothesis, reduced exposure to microbes in early life limits immune system stimulation, potentially increasing the risk of allergic diseases.

End of the analysis — Thank you for joining the workshop!